by William Suzuki

The objective of this work is to build GIFs showing night time luminosity change in some large cities of Brazil. Have a look at the maps below.

Individual Maps

Preliminaries

We start by calling the packages:

library(rgdal) #for readOGR() function
library(raster) #for raster() function
library(tidyverse) #for pipe

The first step to reproduce the maps that we present here is to download the files containing luminosity data: https://ngdc.noaa.gov/eog/dmsp/downloadV4composites.html

For now, I suggest downloading a file for just one year, later we’ll need to download the files for all years. The files are from night time luminosity gathered for each year since 1993 to 2013.

path1 below shows where I stored the files containing raster of luminosity.

#path to directory with raster luminosity night time
path1 = "/home/william/Dropbox/working/RAW_DATA/NightTimeLights/"

Here we create path_luminosity where we’ll store the images generated.

path_luminosity = "/home/william/Dropbox/working/Projects/191130 luminosity Brasil/luminosity_images/"

path_file combines path1 and the place where luminosity data for 1993 is stored.

path_file = paste0(path1,"F101993.v4/F101993.v4b_web.stable_lights.avg_vis.tif")
l1993=raster(path_file) #luminosity 1993

Below we call the Brazillian municipality borders spatial polygons data.

#municipality borders ---------------------------------------------------------
#munic path data 
path_munic = "~/Dropbox/working/RAW_DATA/shpbrasil/municip07.shp"
munic = readOGR(path_munic, layer="municip07")
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/william/Dropbox/working/RAW_DATA/shpbrasil/municip07.shp", layer: "municip07"
## with 5561 features
## It has 3 fields

We can download files for municipality borders from my Dropbox:

https://www.dropbox.com/sh/3ovdzk8erk43aw4/AAAqv4BgWvmZFFWAUIGa7ro2a?dl=0

Note that for the shape file to be read by readOGR() we need to download all files in the folder.

Plot some maps

In this part we explore the sources of data, and plot each map individually the following is an example.

Below is an example for São Paulo city.

The following map is the night time luminosity around the city of São Paulo.

#maps for Sao Paulo 
#long-west, long-east, lat-south, lat-north
bb = c( -47, -46, -24.35,-23.35)
rl1993 = crop(l1993, bb)
plot(rl1993)

This map shows the municipal borders.

munic_crop = crop(munic, bb)
plot(munic_crop)

And this map shows the luminosity and municipality borders together.

plot(rl1993)
lines(munic_crop)

Let’s do just one more example now for the city of Rio de Janeiro.

Note that the object bb below delimits the coordinates for our crop() function. The numbers delimiting it are extremely important I tried to make 1 degree for longitute and latitute. If we want to take the image for another place we’ll need to change bb by hand.

#maps for Rio de Janeiro
#long-west, long-east, lat-south, lat-north
bb = c(-43.8, -42.8, -23.33,-22.33)
rl1993 = crop(l1993, bb)
plot(rl1993)

munic_crop = crop(munic, bb)
plot(munic_crop)

plot(rl1993)
lines(munic_crop)

GIFs

Now we’ll build GIFs for year 1992 to 2013 for many cities of Brazil.

The cities that we analyzed are: Belém, Campo Grande, Fortaleza, Porto Alegre, Salvador, Belo Horizonte, Cuiabá, Goiânia, Recife, São Luis, Brasília, Curitiba, Maceió, Ribeirão Preto, São Paulo, Campinas, Florianopolis, Manaus and Rio de Janeiro.

The package magick we’ll build the GIFs for us.

library(magick) # this is call to animate

The following objects will be important for when we build the loop that will make maps and GIFs for all 22 years and 19 cities.

#path to directory with raster luminosity night time
path1 = "/home/william/Dropbox/working/RAW_DATA/NightTimeLights/"

path_luminosity is the place where we’ll store all *.jpg generated.

path_luminosity = "/home/william/Dropbox/working/Projects/191130 luminosity Brasil/luminosity_images/"

year_vector will be used in the loop.

year_vector = 1992:2013

paths_vector are places where all night time luminosity for each year files are stored.

paths_vector = c(
  "F101992.v4/F101992.v4b_web.stable_lights.avg_vis.tif",
  "F101993.v4/F101993.v4b_web.stable_lights.avg_vis.tif",
  "F101994.v4/F101994.v4b_web.stable_lights.avg_vis.tif",
  "F121995.v4/F121995.v4b_web.stable_lights.avg_vis.tif",
  "F121996.v4/F121996.v4b_web.stable_lights.avg_vis.tif",
  "F141997.v4/F141997.v4b_web.stable_lights.avg_vis.tif",
  "F141998.v4/F141998.v4b_web.stable_lights.avg_vis.tif",
  "F141999.v4/F141999.v4b_web.stable_lights.avg_vis.tif",
  "F152000.v4/F152000.v4b_web.stable_lights.avg_vis.tif",
  "F152001.v4/F152001.v4b_web.stable_lights.avg_vis.tif",
  "F152002.v4/F152002.v4b_web.stable_lights.avg_vis.tif",
  "F152003.v4/F152003.v4b_web.stable_lights.avg_vis.tif",
  "F162004.v4/F162004.v4b_web.stable_lights.avg_vis.tif",
  "F162005.v4/F162005.v4b_web.stable_lights.avg_vis.tif",
  "F162006.v4/F162006.v4b_web.stable_lights.avg_vis.tif",
  "F162007.v4/F162007.v4b_web.stable_lights.avg_vis.tif",
  "F162008.v4/F162008.v4b_web.stable_lights.avg_vis.tif",
  "F162009.v4/F162009.v4b_web.stable_lights.avg_vis.tif",
  "F182010.v4/F182010.v4d_web.stable_lights.avg_vis.tif",
  "F182011.v4/F182011.v4c_web.stable_lights.avg_vis.tif",
  "F182012.v4/F182012.v4c_web.stable_lights.avg_vis.tif",
  "F182013.v4/F182013.v4c_web.stable_lights.avg_vis.tif")

directory_city are the names of the folders where we’ll store the generated *.jpg.

directory_city = 
c("belem","campo_grande",
"fortaleza","porto_alegre","salvador",
"belo_horizonte","cuiaba","goiania",
"recife","sao_luis","brasilia",
"curitiba","maceio","ribeirao_preto",
"sao_paulo","campinas","florianopolis",
"manaus","rio_de_janeiro")

list_bb show the bounding boxes ordered according to the objects above.

list_bb=
list(c(-48.8,-47.8,-1.8,-0.8),
c(-55,-54,-20.9,-19.9),
c(-39.1,-38.1,-4.3,-3.3),
c(-51.5,-50.5,-30.4,-29.4),
c(-39.2,-38.2,-13.2,-12.2),
c(-44.5,-43.5,-20.3,-19.3),
c(-56.7,-55.7,-16,-15),
c(-49.7,-48.7,-17.1,-16.1),
c(-35.7,-34.7, -8.6,-7.6),
c(-44.7,-43.7,-3.2,-2.2),
c(-48.3,-47.3,-16.3,-15.3),
c(-49.7,-48.7,-25.9,-24.9),
c(-36.2,-35.2,-10,-9),
c(-48.2,-47.2,-21.5,-20.5),
c( -47, -46, -24.35,-23.35),
c(-47.5,-46.5,-23.2,-22.2),
c(-49,-48,-28,-27),
c(-60.4,-59.4,-3.5,-2.5),
c(-43.8, -42.8, -23.33,-22.33))

Now we build the loop. This loop is building for each year, for each city a .jpg file and storing it in the appropriate folder. The in the end the loop builds a GIF for each city.

#loop ---------------------------------------------------------
#1:19 all cities
jj = 1
for (jj in 1:19) { #loop for city
  bb = list_bb[[jj]]
  munic_crop = crop(munic, bb)
  
  ii = 1
  for (ii in 1:22){ #loop for year
    path_file = paste0(path1,paths_vector[ii])
    lumiRaster=raster(path_file)
    lum_crop = crop(lumiRaster, bb)
    
    xlab_title = paste(directory_city[jj], year_vector[ii])
    
    path_save = paste0(path_luminosity,directory_city[jj],"/","lum",as.character(year_vector[ii]),".jpg")
    jpeg(path_save)
    plot(lum_crop, xlab=xlab_title)
    lines(munic_crop)
    dev.off() 
    print(xlab_title)
  }
  path_city = paste0(path_luminosity,directory_city[jj])
  list.files(path = path_city, pattern = "*.jpg", full.names = T) %>% 
    map(image_read) %>% # reads each path file
    image_join() %>% # joins image
    image_animate(fps=1) %>% # animates, can opt for number of loops
    image_write(paste0(directory_city[jj],"Lum.gif")) # write to current dir
}

Here we show the GIFs generated: